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Abstract . The determination of very precise 
orbits and geodynamic parameters from laser 
tracking data requires the continual development 
and improvement of the software systems and 
computational techniques . Computational accura- 
cies at the few centimeter level are presently 
required to match the performance of the present 
day laser ranging systems and altimeters and in 
the next few years the accuracies are expected 
to increase further. In this paper and the major 
error sources in orbit determination are briefly 
discussed and the present and future modeling 
activities needed to meet the accuracy require- 
ments of the next few years are described. 

Introduction 

The precise computation of the motion of earth 
satellites has become a critical component of 
several space techniques for studying the earth . 
Perhaps the two most relevant at the present time 
are the computation of orbits for the determina- 
tion of crustal motions and the computation of 
orbits for altimeter satellites . These and other 
applications require orbital accuracies in the 
centimeter range and challenge our present-day 
capability to accurately model many of the forces 
that perturb satellites and also our computational 
techniques . Laser ranging systems operated by 
NASA Goddard Space Flight Center are now at about 
the 5 cm accuracy level and systems capable of 
2 cm are presently under development (Silverberg, 
this volume) . However, our present ability to 
determine the orbits of satellites has not yet 
achieved the same level of accuracy and histori- 
cally orbital accuracy has lagged behind the 
observational accuracy by several years because 
the improved observations are needed to improve 
upon the models for the spacecraft behavior. 

This situation is particularly true for the 
gravity field, which for most satellites and 
applications is still the largest source of 
orbital error. 

Figure 1 shows the developmetn of orbit deter- 
mination for arc lengths of three to five revolu- 
tions over the last several years in comparison 
with the improvement in the quality of laser 
range measurements . Figure 1 represents the 
situation for a typical low altitude spacecraft; 
such as Beacon Explorer C or GEOS -3 which have 
altitudes of about 1000 km or less . For most 
higher altitude satellites, such as Lageos at 
6000 km, the orbital fit will be closer to the 
data quality. Another important consideration 
in interpreting Figure 1 is that the orbit fit 
curve is largely based on the experience in the 
San Andreas fault experiment (Smith, et al. f 
1979) in which two or three laser tracking sta^ 
tions within one or two thousand kilometers of 
each other were used. If these stations had 
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been on different continents the fit might have 
been larger . 

For these short orbital arcs of a few revolu- 
tions (3 to 8 hours) the largest source of error 
in the orbital computations is usually the gravity 
field but as the arc length increases to days , 
weeks, and months other sources of error usually 
begin to dominate. Figure 2 illustrates the way 
the orbit error typically grows as a function of 
arc length for the major perturbing forces. The 
vertical axis is arbitrary in scale and is only 
representative of the relative magnitudes. The 
exact variation of each of the curves will depend 
on the satellite, its orbit and the sophistica- 
tion of the modeling of the perturbing force used 
in the orbit computation. For example, numerical 
errors in the computations are usually negligible 
if an analytical theory is being used for the 
spacecraft motion but can become extremely impor- 
tant if a numerical integration system is used 
because the error accumulates as the arc length 
increases . 

At GSFC the major software system for the 
computation of precision orbits is the Geodyn 
program (Putney, 1977) . This program system has 
the ability to determine orbits from a variety 
of tracking data types and is capable of esti- 
mating various geophysical parameters such as 
polar motion and earth rotation, tidal para- 
meters, geopotential coefficients, as well as 
parameters associated with the perturbing forces , 
such as drag coefficients, and measurement errors, 
such as range biases. 

In the following section the status of our 
modeling of the different perturbing forces at 
GSFC will be briefly described together with our 
plans for future improvements . 

Orbital Perturbing Forces 

Gravity 

One of the areas that has shown greatest 
improvement during the last decade has been our 
knowledge of the gravity field of the earth. The 
inclusion of high precision range measurements 
from laser tracking systems and Unified S-Band 
doppler data and, more recently, altimeter data 
has enabled our definition of the gravity field 
to extend out to degree and order 36 with 
specialized altimeter and surface gravity solu- 
tions extending out to degree and order 180 . 

These new fields (Goddard Earth Models - GEM) 
developed at GSFC have permitted improvements of 
at least an order of magnitude in the determina^ 
tion of short orbital arcs of satellites over the 
past decade . Figure 3 shows a comparison of the 
abilities of three gravity models, GEM's 1, 7 
and 9, to fit five consecutive passes of laser 
data from a single tracking station. These five 
passes, obtained at GSFC in 1974 on the Beacon 
Explorer C satellite, when analyzed by the GEM 1 
gravity field developed in 1970-71 (Smith, et al. , 
1973) could only be satisfied at the 2 meter 
level even though the data was of 10 cm quality . 
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The same data analyzed a few years later with 
ghe GEM 7 model (Wagner, et al., 1977) could be 
satisfied to about the 50 centimeter level and 
more recently the GEM 9 (Lerch, et al., 1978) 
model fits to 12 cm. The improvement from GEM 1 
to GEM 9 has been brought about largely by the 
inclusion in the later models of large quantities 
of laser tracking data on several satellites, 
particularly GEOS-3, but not the Beacon Explorer 
C data shown in Figure 3 . The slight curvature 
of the GEM 9 results in Figure 3 show that some 
gravitational signature still exists in the data 
and that some improvement still remains to be 
made although this may well be the most difficult. 

For longer arcs the gravity error increases 
to about the 50 cm to 1 meter level after one 
week with the GEM 9 field for low altitude satel- 
lites such as Lageos, reaching 50 cm after about 
one month. 

The present plans at GSFC are to continue to 
improve our knowledge of the gravity field so 
that the locations of mobile and fixed laser sta- 
tions can be determined to the few centimeter 
level for the precise measurement of crustal 
motions; and also for the precise analysis of the 
GEOS-3 and Seasat altimeter data. 

Air Drag 

As the orbital arc length increases the second 
most important perturbing force (after gravity) 
for low altitude satellites is usually the effect 
of air drag. The general form of the perturba- 
tion is 


porated by computing externally to the main pro- 
gram the cross-sectional area as a function of 
angle of incidence and including this information 
in tabular form in the orbit program. At each 
integration step the appropriate area is deduced 
from the table and used in the drag calculations. 

Solar Radiation Pressure 

The form of the perturbation by solar radia- 
tion pressure is 

_ A (solar constant) 


where C R is a constant that accounts for the type 
of reflection, specular diffuse, etc., that is 
taking place at the satellite, ^ is the area (A) 
to mass (M) ratio and C is the velocity of light. 
The model used in the Geodyn program includes a 
solar flux varying with solar distance, and 
approximations for absorption and refraction at 
the terminator . The incident area is variable 
in the same manner as for air drag and takes into 
account spacecraft attitude, shape, shadowing 
and varying reflective properties over the 
spacecraft. Although these computations are 
reasonably precise we believe that errors are 
occurring at the umbra/penumbra/full sunlight 
boundaries where the numerical integration 
procedure jumps over one or even two boundaries 
in one step. The effect of this error is esti- 
mated to be a slight offset in the boundary 
location . 


accel . 



Di P 


V 


2 


where C p is the drag coefficient, A is the space- 
craft cross-sectional area, M is the spacecraft 
mass, p is the air density and v is the space- 
craft velocity. At the present time the density 
models used in computing the drag acceleration 
are based on the work of Jacchia and include 
variations in solar activity, diurnal terms, 
geomagnetic effects and semi-annual and seasonal 
latitude variations. In order to improve the 
responsiveness of the model to unmodeled changes 
in density we have introduced a time dependent 
parameter (p) that enables us to account for 
systemated changes during the orbital arc . In 
addition, we are introducing a capability to 
estimate the drag coefficients (C ) for specific 
periods during the orbital arc. Thus it will be 
possible to vary the drag coefficient from one 
day to the next and thereby modify the drag ac- 
celeration without any change to the density 
model. We believe this may accomodate density 
changes that last for short periods of time that 
are not represented in our models . 

Another capability that exists in the GSFC 
Geodyn program is to accurately model the cross- 
sectional area of the spacecraft. This facility 
was introduced for the Beacon Explorer C, GEOS-3 
and Seasat spacecraft because of their irregu- 
larity in shape and the need for very precise 
orbit calculations on these satellites. In all 
these cases drag was a major influence in the 
orbit computations and inclusion of a variable 
cross-sectional area could improve upon the 
computational accuracy. The technique is incor- 


Albedo Radiation Pressure 

At the present time we do not have an albedo 
model in our orbit computation program but one 
is under consideration. The basic form of the 
perturbation is 
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and is similar to that for solar radiation pres- 
sure except for the (R/D) term which shows that 
the acceleration follows the inverse square law 
(R is the earth radius, D is the radial distance 
of the satellite) . The difficulty with evalua- 
ting this perturbing force is that the albedo is 
variable in both space and time and needs to be 
evaluated at every integration step for the en- 
tire surface of the earth observable from the 
spacecraft. This procedure is computationally 
very time consuming. Simplification of the 
albedo into only day/night effects , for example , 
will probably underestimate the effect and pro- 
vide deceiving results . Our present considera- 
tions are directed towards the computation of 
the long-term effects of albedo by digitally 
representing seasonal albedo maps of the world 
derived from satellite meteorological data. 


Earth and Ocean Tides 

Our present modeling of earth tides in the 
Geodyn program is a single second degree spheri- 
cal harmonic with one amplitude and phase. This 
model is used to compute the gravitational effects 
of earth tides on the satellite orbit and to 
compute the body tides on the locations of the 
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tracking stations. Our ocean tide models only 
account for the displacement of the ocean surface 
(Hendershott model globally, Mofjeld in Northwest 
Atlantic) for the analysis of satellite altimeter 
data and do not include the gravitational effect 
on the satellite. However, considerable accomo- 
dation of the ocean tide effect on the orbit can 
be achieved by suitable modification of the solid 
earth tidal amplitude and phase . For example , 
for Beacon Explorer C we were able to model 
approximately 90 to 95% of the combined long 
period earth and ocean tidal perturbation of the 
satellite with a Love number (k 2 ) of 0.245, and 
a phase lag of 3.2 degrees (Smith, et al., 1972) 
used in the solid earth tide model. This 
accomodation of the oceans within the earth model 
ignores any frequency dependent terms in either 
the earth or ocean tides. Consequently, we are 
planning to allow for a frequency dependence of 
both the amplitude and phase of the solid-earth 
which we expect to permit complete accomodation 
of the ocean tides. In addition, we expect to 
incorporate a spherical harmonic representation 
of the major ocean tides (M 2 , S^, K , K,, , N^, 

O^, P^, etc), derived from the numerical integra- 
tion of Laplace's tidal equations, into our 
program system so that we can use these models 
for both altimetry and orbital analysis and also 
be able to use altimetry and orbital data to 
improve on the coefficients in the tidal 
expansions. 

Numerical Problems 

Numerical integration systems introduce errors 
of rounding and truncation into the orbit calcu- 
lations that can become significant for long 
orbital arcs. In the Geodyn program system a 
typical step size within the integrator will be 
about sixty seconds but as the size of the gra- 
vity models has increased we are finding that 
this figure needs to be reduced to perhaps forty 
seconds in order to properly account for the high 
frequency terms; and the CPU time increases 
accordingly. Generally, with step sizes of the 
order of sixty seconds integration error can be 
kept to the order of a meter after about 30 days . 
This error is also predominantly along track and 
is an acceleration similar to drag. Thus, if 
drag is being adjusted in the orbit determination 
process then the integration error will be largely 
absorbed into the adjustment. In this case inte- 
gration error is not a major contributor to orbit 
error. 

Another aspect of numerical problems in orbit 
determination is the core size required to deter- 
mine large spherical harmonic gravity fields and 
associated station coordinates. Some of our 
large gravity field and station coordinate solu- 
tions have nearly 2,000 unknowns and require 
several million bytes on an IBM 360/95 type 
computer. This storage is not always available 
and if the computations are attempted in smaller 
core the CPU and/or 10 time increases . Nearly 
all our computations in satellite orbit, geodesy 
and geodynamics are requiring greater accuracy 
today than a few years ago and this means an 
increase in both CPU time and core storage. 

Because of this need for additional precision 
and complexity in the computations we are giving 
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. Fig. 1. Development of orbit determination 
capability . 

consideration to the advantages and disadvantages 
of parallel processor computers and dedicated 
minicomputers for some of this work . 

Conclusions 

The application of satellite geodetic techni- 
ques to problems in solid earth and ocean 
dynamics is requiring ever increasing accuracy in 
the computation of satellite orbits . This need 
presents considerable difficulty in the modeling 
of the many perturbing forces that influence the 
spacecraft motion. With increasing complexity 
and accuracy is the need for faster and larger 
computing facilities. 

In this paper I have briefly described the 
status of the major GSFC orbit and geodynamic 
parameter estimation program (Geodyn) and the 
degree of complexity that we are finding neces- 
sary to meet orbital and geophysical accuracies. 
The computing of precision orbits at the centi- 
meter level is difficult to obtain and even more 
difficult to maintain for any length of time and 
may have to be limited to satellites, such as 
Lageos, that are carefully configured to mini- < 
mize the perturbing forces. 
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